#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstring>
#include <string>
#include <iostream>
#include <fstream>
#include <unistd.h>
#include <sys/socket.h>
#include <sys/types.h>
#include <netinet/in.h>
#include <netdb.h>
#include <sys/time.h>
#include <sys/types.h>
#include <sys/sysinfo.h>
#include <cstdlib>
#include <cfloat>
#include <cstddef>
#include <stdexcept>
#include <cstdio>
#include <map>
#include <sstream>
#include <cctype>
#include <iomanip>
#include <sstream>
#include <vector>
#include <limits>

// jsonxx versioning: major.minor-extra where
// major = { number }
// minor = { number }
// extra = { 'a':alpha, 'b':beta, 'rc': release candidate, 'r': release, 's':stable }
#define JSONXX_MAJOR    "0"
#define JSONXX_MINOR    "22"
#define JSONXX_EXTRA    "a"
#define JSONXX_VERSION  JSONXX_MAJOR "." JSONXX_MINOR "-" JSONXX_EXTRA
#define JSONXX_XML_TAG  "<!-- generated by jsonxx " JSONXX_VERSION " -->"

#if __cplusplus > 199711L
#define JSONXX_COMPILER_HAS_CXX11 1
#elif defined(_MSC_VER) && _MSC_VER > 1700
#define JSONXX_COMPILER_HAS_CXX11 1
#else
#define JSONXX_COMPILER_HAS_CXX11 0
#endif

#define JSONXX_ASSERT(...) do { if( jsonxx::Assertions ) \
  jsonxx::assertion(__FILE__,__LINE__,#__VA_ARGS__,bool(__VA_ARGS__)); } while(0)

namespace jsonxx {

// Settings
enum Settings {
  // constants
  Enabled = true,
  Disabled = false,
  Permissive = true,
  Strict = false,
  // values
  Parser = Permissive,  // permissive or strict parsing
  UnquotedKeys = Disabled, // support of unquoted keys
  Assertions = Enabled  // enabled or disabled assertions (these asserts work both in DEBUG and RELEASE builds)
};

enum Format {
  JSON      = 0,     // JSON output
  JSONx     = 1,     // XML output, JSONx  format. see http://goo.gl/I3cxs
};

// Types
typedef long double Number;
typedef bool Boolean;
typedef std::string String;
struct Null {};
class Value;
class Object;
class Array;

// Identity meta-function
template <typename T>
struct identity {
  typedef T type;
};

// Detail
void assertion( const char *file, int line, const char *expression, bool result );

// A JSON Object
class Object {
 public:
  Object();
  ~Object();

  template <typename T>
  bool has(const std::string& key) const;

  // Always call has<>() first. If the key doesn't exist, consider
  // the behavior undefined.
  template <typename T>
  T& get(const std::string& key);
  template <typename T>
  const T& get(const std::string& key) const;

  template <typename T>
  const T& get(const std::string& key, const typename identity<T>::type& default_value) const;

  size_t size() const;
  bool empty() const;

  const std::map<std::string, Value*>& kv_map() const;

  void reset();
  bool parse(std::istream &input);
  bool parse(const std::string &input);
  typedef std::map<std::string, Value*> container;
  void import( const Object &other );
  void import( const std::string &key, const Value &value );
  Object &operator<<(const Value &value);
  Object &operator<<(const Object &value);
  Object &operator=(const Object &value);
  Object(const Object &other);
  Object(const std::string &key, const Value &value);
  template<size_t N>
  Object(const char (&key)[N], const Value &value) {
    import(key,value);
  }
  template<typename T>
  Object &operator<<(const T &value);

 protected:
  static bool parse(std::istream& input, Object& object);
  container value_map_;
  std::string odd;
};

class Array {
 public:
  Array();
  ~Array();

  size_t size() const;
  bool empty() const;

  template <typename T>
  bool has(unsigned int i) const;

  template <typename T>
  T& get(unsigned int i);
  template <typename T>
  const T& get(unsigned int i) const;

  template <typename T>
  const T& get(unsigned int i, const typename identity<T>::type& default_value) const;

  const std::vector<Value*>& values() const {
    return values_;
  }

  void reset();
  bool parse(std::istream &input);
  bool parse(const std::string &input);
  typedef std::vector<Value*> container;
  void import(const Array &other);
  void import(const Value &value);
  Array &operator<<(const Array &other);
  Array &operator<<(const Value &value);
  Array &operator=(const Array &other);
  Array &operator=(const Value &value);
  Array(const Array &other);
  Array(const Value &value);
 protected:
  static bool parse(std::istream& input, Array& array);
  container values_;
};

// A value could be a number, an array, a string, an object, a
// boolean, or null
class Value {
 public:

  Value();
  ~Value() { reset(); }
  void reset();

  template<typename T>
  void import( const T & ) {
    reset();
    type_ = INVALID_;
    // debug
    // std::cout << "[WARN] No support for " << typeid(t).name() << std::endl;
  }
  void import( const bool &b ) {
    reset();
    type_ = BOOL_;
    bool_value_ = b;
  }
#define $number(TYPE) \
  void import( const TYPE &n ) { \
    reset(); \
    type_ = NUMBER_; \
    number_value_ = static_cast<long double>(n); \
  }
  $number( char )
  $number( int )
  $number( long )
  $number( long long )
  $number( unsigned char )
  $number( unsigned int )
  $number( unsigned long )
  $number( unsigned long long )
  $number( float )
  $number( double )
  $number( long double )
#undef $number
#if JSONXX_COMPILER_HAS_CXX11 > 0
  void import( const std::nullptr_t & ) {
    reset();
    type_ = NULL_;
  }
#endif
  void import( const Null & ) {
    reset();
    type_ = NULL_;
  }
  void import( const String &s ) {
    reset();
    type_ = STRING_;
    *( string_value_ = new String() ) = s;
  }
  void import( const Array &a ) {
    reset();
    type_ = ARRAY_;
    *( array_value_ = new Array() ) = a;
  }
  void import( const Object &o ) {
    reset();
    type_ = OBJECT_;
    *( object_value_ = new Object() ) = o;
  }
  void import( const Value &other ) {
    if (this != &other)
    switch (other.type_) {
      case NULL_:
        import( Null() );
        break;
      case BOOL_:
        import( other.bool_value_ );
        break;
      case NUMBER_:
        import( other.number_value_ );
        break;
      case STRING_:
        import( *other.string_value_ );
        break;
      case ARRAY_:
        import( *other.array_value_ );
        break;
      case OBJECT_:
        import( *other.object_value_ );
        break;
      case INVALID_:
        type_ = INVALID_;
        break;
      default:
        JSONXX_ASSERT( !"not implemented" );
    }
  }
  template<typename T>
  Value &operator <<( const T &t ) {
    import(t);
    return *this;
  }
  template<typename T>
  Value &operator =( const T &t ) {
    reset();
    import(t);
    return *this;
  }
  Value(const Value &other);
  template<typename T>
  Value( const T&t ) : type_(INVALID_) { import(t); }
  template<size_t N>
  Value( const char (&t)[N] ) : type_(INVALID_) { import( std::string(t) ); }

  bool parse(std::istream &input);
  bool parse(const std::string &input);

  template<typename T>
  bool is() const;
  template<typename T>
  T& get();
  template<typename T>
  const T& get() const;

  bool empty() const;

 public:
  enum {
    NUMBER_,
    STRING_,
    BOOL_,
    NULL_,
    ARRAY_,
    OBJECT_,
    INVALID_
  } type_;
  union {
    Number number_value_;
    String* string_value_;
    Boolean bool_value_;
    Array* array_value_;
    Object* object_value_;
  };

protected:
  static bool parse(std::istream& input, Value& value);
};

template <typename T>
bool Array::has(unsigned int i) const {
  if (i >= size()) {
    return false;
  } else {
    Value* v = values_.at(i);
    return v->is<T>();
  }
}

template <typename T>
T& Array::get(unsigned int i) {
  JSONXX_ASSERT(i < size());
  Value* v = values_.at(i);
  return v->get<T>();
}

template <typename T>
const T& Array::get(unsigned int i) const {
  JSONXX_ASSERT(i < size());
  const Value* v = values_.at(i);
  return v->get<T>();
}

template <typename T>
const T& Array::get(unsigned int i, const typename identity<T>::type& default_value) const {
  if(has<T>(i)) {
    const Value* v = values_.at(i);
    return v->get<T>();
  } else {
    return default_value;
  }
}

template <typename T>
bool Object::has(const std::string& key) const {
  container::const_iterator it(value_map_.find(key));
  return it != value_map_.end() && it->second->is<T>();
}

template <typename T>
T& Object::get(const std::string& key) {
  JSONXX_ASSERT(has<T>(key));
  return value_map_.find(key)->second->get<T>();
}

template <typename T>
const T& Object::get(const std::string& key) const {
  JSONXX_ASSERT(has<T>(key));
  return value_map_.find(key)->second->get<T>();
}

template <typename T>
const T& Object::get(const std::string& key, const typename identity<T>::type& default_value) const {
  if (has<T>(key)) {
    return value_map_.find(key)->second->get<T>();
  } else {
    return default_value;
  }
}
    
template<>
inline bool Value::is<Value>() const {
    return true;
}

template<>
inline bool Value::is<Null>() const {
  return type_ == NULL_;
}

template<>
inline bool Value::is<Boolean>() const {
  return type_ == BOOL_;
}

template<>
inline bool Value::is<String>() const {
  return type_ == STRING_;
}

template<>
inline bool Value::is<Number>() const {
  return type_ == NUMBER_;
}

template<>
inline bool Value::is<Array>() const {
  return type_ == ARRAY_;
}

template<>
inline bool Value::is<Object>() const {
  return type_ == OBJECT_;
}
    
template<>
inline Value& Value::get<Value>() {
    return *this;
}
    
template<>
inline const Value& Value::get<Value>() const {
    return *this;
}

template<>
inline bool& Value::get<Boolean>() {
  JSONXX_ASSERT(is<Boolean>());
  return bool_value_;
}

template<>
inline std::string& Value::get<String>() {
  JSONXX_ASSERT(is<String>());
  return *string_value_;
}

template<>
inline Number& Value::get<Number>() {
  JSONXX_ASSERT(is<Number>());
  return number_value_;
}

template<>
inline Array& Value::get<Array>() {
  JSONXX_ASSERT(is<Array>());
  return *array_value_;
}

template<>
inline Object& Value::get<Object>() {
  JSONXX_ASSERT(is<Object>());
  return *object_value_;
}

template<>
inline const Boolean& Value::get<Boolean>() const {
  JSONXX_ASSERT(is<Boolean>());
  return bool_value_;
}

template<>
inline const String& Value::get<String>() const {
  JSONXX_ASSERT(is<String>());
  return *string_value_;
}

template<>
inline const Number& Value::get<Number>() const {
  JSONXX_ASSERT(is<Number>());
  return number_value_;
}

template<>
inline const Array& Value::get<Array>() const {
  JSONXX_ASSERT(is<Array>());
  return *array_value_;
}

template<>
inline const Object& Value::get<Object>() const {
  JSONXX_ASSERT(is<Object>());
  return *object_value_;
}

template<typename T>
inline Object &Object::operator<<(const T &value) {
  return *this << Value(value), *this;
}

}  // namespace jsonxx

#if defined(NDEBUG) || defined(_NDEBUG)
#   define JSONXX_REENABLE_NDEBUG
#   undef  NDEBUG
#   undef _NDEBUG
#endif
void jsonxx::assertion( const char *file, int line, const char *expression, bool result ) {
    if( !result ) {
        fprintf( stderr, "[JSONXX] expression '%s' failed at %s:%d -> ", expression, file, line );
        assert( 0 );
    }
}
#if defined(JSONXX_REENABLE_NDEBUG)
#   define  NDEBUG
#   define _NDEBUG
#endif
namespace jsonxx {

//static_assert( sizeof(unsigned long long) < sizeof(long double), "'long double' cannot hold 64bit values in this compiler :(");

bool match(const char* pattern, std::istream& input);
bool parse_array(std::istream& input, Array& array);
bool parse_bool(std::istream& input, Boolean& value);
bool parse_comment(std::istream &input);
bool parse_null(std::istream& input);
bool parse_number(std::istream& input, Number& value);
bool parse_object(std::istream& input, Object& object);
bool parse_string(std::istream& input, String& value);
bool parse_identifier(std::istream& input, String& value);
bool parse_value(std::istream& input, Value& value);

// Try to consume characters from the input stream and match the
// pattern string.
bool match(const char* pattern, std::istream& input) {
    input >> std::ws;
    const char* cur(pattern);
    char ch(0);
    while(input && !input.eof() && *cur != 0) {
        input.get(ch);
        if (ch != *cur) {
            input.putback(ch);
            if( parse_comment(input) )
                continue;
            while (cur > pattern) {
                cur--;
                input.putback(*cur);
            }
            return false;
        } else {
            cur++;
        }
    }
    return *cur == 0;
}

bool parse_string(std::istream& input, String& value) {
    char ch = '\0', delimiter = '"';
    if (!match("\"", input))  {
        if (Parser == Strict) {
            return false;
        }
        delimiter = '\'';
        if (input.peek() != delimiter) {
            return false;
        }
        input.get(ch);
    }
    while(!input.eof() && input.good()) {
        input.get(ch);
        if (ch == delimiter) {
            break;
        }
        if (ch == '\\') {
            input.get(ch);
            switch(ch) {
                case '\\':
                case '/':
                    value.push_back(ch);
                    break;
                case 'b':
                    value.push_back('\b');
                    break;
                case 'f':
                    value.push_back('\f');
                    break;
                case 'n':
                    value.push_back('\n');
                    break;
                case 'r':
                    value.push_back('\r');
                    break;
                case 't':
                    value.push_back('\t');
                    break;
                case 'u': {
                        int i;
                        std::stringstream ss;
                        for( i = 0; (!input.eof() && input.good()) && i < 4; ++i ) {
                            input.get(ch);
                            ss << std::hex << ch;
                        }
                        if( input.good() && (ss >> i) )
                            value.push_back(i);
                    }
                    break;
                default:
                    if (ch != delimiter) {
                        value.push_back('\\');
                        value.push_back(ch);
                    } else value.push_back(ch);
                    break;
            }
        } else {
            value.push_back(ch);
        }
    }
    if (input && ch == delimiter) {
        return true;
    } else {
        return false;
    }
}

bool parse_identifier(std::istream& input, String& value) {
    input >> std::ws;

    char ch = '\0', delimiter = ':';
    bool first = true;

    while(!input.eof() && input.good()) {
        input.get(ch);

        if (ch == delimiter) {
            input.unget();
            break;
        }

        if(first) {
            if ((ch != '_' && ch != '$') &&
                    (ch < 'a' || ch > 'z') &&
                    (ch < 'A' || ch > 'Z')) {
                return false;
            }
            first = false;
        }
        if(ch == '_' || ch == '$' ||
            (ch >= 'a' && ch <= 'z') ||
            (ch >= 'A' && ch <= 'Z') ||
            (ch >= '0' && ch <= '9')) {
            value.push_back(ch);            
        }
        else if(ch == '\t' || ch == ' ') {
            input >> std::ws;
        }
    }
    if (input && ch == delimiter) {
        return true;
    } else {
        return false;
    }
}

bool parse_number(std::istream& input, Number& value) {
    input >> std::ws;
    std::streampos rollback = input.tellg();
    input >> value;
    if (input.fail()) {
        input.clear();
        input.seekg(rollback);
        return false;
    }
    return true;
}

bool parse_bool(std::istream& input, Boolean& value) {
    if (match("true", input))  {
        value = true;
        return true;
    }
    if (match("false", input)) {
        value = false;
        return true;
    }
    return false;
}

bool parse_null(std::istream& input) {
    if (match("null", input))  {
        return true;
    }
    if (Parser == Strict) {
        return false;
    }
    return (input.peek()==',');
}

bool parse_array(std::istream& input, Array& array) {
    return array.parse(input);
}

bool parse_object(std::istream& input, Object& object) {
    return object.parse(input);
}

bool parse_comment(std::istream &input) {
    if( Parser == Permissive )
    if( !input.eof() && input.peek() == '/' )
    {
        char ch0(0);
        input.get(ch0);

        if( !input.eof() )
        {
            char ch1(0);
            input.get(ch1);

            if( ch0 == '/' && ch1 == '/' )
            {
                // trim chars till \r or \n
                for( char ch(0); !input.eof() && (input.peek() != '\r' && input.peek() != '\n'); )
                    input.get(ch);

                // consume spaces, tabs, \r or \n, in case no eof is found
                if( !input.eof() )
                    input >> std::ws;
                return true;
            }

            input.unget();
            input.clear();
        }

        input.unget();
        input.clear();
    }

    return false;
}

bool parse_value(std::istream& input, Value& value) {
    return value.parse(input);
}


Object::Object() : value_map_() {}

Object::~Object() {
    reset();
}

bool Object::parse(std::istream& input, Object& object) {
    object.reset();

    if (!match("{", input)) {
        return false;
    }
    if (match("}", input)) {
        return true;
    }

    do {
        std::string key;
        if(UnquotedKeys == Enabled) {
            if (!parse_identifier(input, key)) {
                if (Parser == Permissive) {
                    if (input.peek() == '}')
                        break;
                }
                return false;
            }
        }
        else {
            if (!parse_string(input, key)) {
                if (Parser == Permissive) {
                    if (input.peek() == '}')
                        break;
                }
                return false;
            }
        }
        if (!match(":", input)) {
            return false;
        }
        Value* v = new Value();
        if (!parse_value(input, *v)) {
            delete v;
            break;
        }
        object.value_map_[key] = v;
    } while (match(",", input));


    if (!match("}", input)) {
        return false;
    }

    return true;
}

Value::Value() : type_(INVALID_) {}

void Value::reset() {
    if (type_ == STRING_) {
        delete string_value_;
        string_value_ = 0;
    }
    else if (type_ == OBJECT_) {
        delete object_value_;
        object_value_ = 0;
    }
    else if (type_ == ARRAY_) {
        delete array_value_;
        array_value_ = 0;
    }
}

bool Value::parse(std::istream& input, Value& value) {
    value.reset();

    std::string string_value;
    if (parse_string(input, string_value)) {
        value.string_value_ = new std::string();
        value.string_value_->swap(string_value);
        value.type_ = STRING_;
        return true;
    }
    if (parse_number(input, value.number_value_)) {
        value.type_ = NUMBER_;
        return true;
    }

    if (parse_bool(input, value.bool_value_)) {
        value.type_ = BOOL_;
        return true;
    }
    if (parse_null(input)) {
        value.type_ = NULL_;
        return true;
    }
    if (input.peek() == '[') {
        value.array_value_ = new Array();
        if (parse_array(input, *value.array_value_)) {
            value.type_ = ARRAY_;
            return true;
        }
        delete value.array_value_;
    }
    value.object_value_ = new Object();
    if (parse_object(input, *value.object_value_)) {
        value.type_ = OBJECT_;
        return true;
    }
    delete value.object_value_;
    return false;
}

Array::Array() : values_() {}

Array::~Array() {
    reset();
}

bool Array::parse(std::istream& input, Array& array) {
    array.reset();

    if (!match("[", input)) {
        return false;
    }
    if (match("]", input)) {
        return true;
    }

    do {
        Value* v = new Value();
        if (!parse_value(input, *v)) {
            delete v;
            break;
        }
        array.values_.push_back(v);
    } while (match(",", input));

    if (!match("]", input)) {
        return false;
    }
    return true;
}

}  // namespace jsonxx


namespace jsonxx {
namespace {

typedef unsigned char byte;

//template<bool quote>
std::string escape_string( const std::string &input, const bool quote = false ) {
    static std::string map[256], *once = 0;
    if( !once ) {
        // base
        for( int i = 0; i < 256; ++i ) {
            map[ i ] = std::string() + char(i);
        }
        // non-printable
        for( int i = 0; i < 32; ++i ) {
            std::stringstream str;
            str << "\\u" << std::hex << std::setw(4) << std::setfill('0') << i;
            map[ i ] = str.str();
        }
        // exceptions
        map[ byte('"') ] = "\\\"";
        map[ byte('\\') ] = "\\\\";
        map[ byte('/') ] = "\\/";
        map[ byte('\b') ] = "\\b";
        map[ byte('\f') ] = "\\f";
        map[ byte('\n') ] = "\\n";
        map[ byte('\r') ] = "\\r";
        map[ byte('\t') ] = "\\t";

        once = map;
    }
    std::string output;
    output.reserve( input.size() * 2 + 2 ); // worst scenario
    if( quote ) output += '"';
    for( std::string::const_iterator it = input.begin(), end = input.end(); it != end; ++it )
        output += map[ byte(*it) ];
    if( quote ) output += '"';
    return output;
}


namespace json {

    std::string remove_last_comma( const std::string &_input ) {
        std::string input( _input );
        size_t size = input.size();
        if( size > 2 )
            if( input[ size - 2 ] == ',' )
                input[ size - 2 ] = ' ';
        return input;
    }

    std::string tag( unsigned format, unsigned depth, const std::string &name, const jsonxx::Value &t) {
        std::stringstream ss;
        const std::string tab(depth, '\t');

        if( !name.empty() )
            ss << tab << '\"' << escape_string( name ) << '\"' << ':' << ' ';
        else
            ss << tab;

        switch( t.type_ )
        {
            default:
            case jsonxx::Value::NULL_:
                ss << "null";
                return ss.str() + ",\n";

            case jsonxx::Value::BOOL_:
                ss << ( t.bool_value_ ? "true" : "false" );
                return ss.str() + ",\n";

            case jsonxx::Value::ARRAY_:
                ss << "[\n";
                for(Array::container::const_iterator it = t.array_value_->values().begin(),
                    end = t.array_value_->values().end(); it != end; ++it )
                  ss << tag( format, depth+1, std::string(), **it );
                return remove_last_comma( ss.str() ) + tab + "]" ",\n";

            case jsonxx::Value::STRING_:
                ss << '\"' << escape_string( *t.string_value_ ) << '\"';
                return ss.str() + ",\n";

            case jsonxx::Value::OBJECT_:
                ss << "{\n";
                for(Object::container::const_iterator it=t.object_value_->kv_map().begin(),
                    end = t.object_value_->kv_map().end(); it != end ; ++it)
                  ss << tag( format, depth+1, it->first, *it->second );
                return remove_last_comma( ss.str() ) + tab + "}" ",\n";

            case jsonxx::Value::NUMBER_:
                // max precision
                ss << std::setprecision(std::numeric_limits<long double>::digits10 + 1);
                ss << t.number_value_;
                return ss.str() + ",\n";
        }
    }
} // namespace jsonxx::anon::json


} // namespace jsonxx::anon

Object::Object(const Object &other) {
  import(other);
}
Object::Object(const std::string &key, const Value &value) {
  import(key,value);
}
void Object::import( const Object &other ) {
  odd.clear();
  if (this != &other) {
    // default
    container::const_iterator
        it = other.value_map_.begin(),
        end = other.value_map_.end();
    for (/**/ ; it != end ; ++it) {
      container::iterator found = value_map_.find(it->first);
      if( found != value_map_.end() ) {
        delete found->second;
      }
      value_map_[ it->first ] = new Value( *it->second );
    }
  } else {
    // recursion is supported here
    import( Object(*this) );
  }
}
void Object::import( const std::string &key, const Value &value ) {
  odd.clear();
  container::iterator found = value_map_.find(key);
  if( found != value_map_.end() ) {
    delete found->second;
  }
  value_map_[ key ] = new Value( value );
}
Object &Object::operator=(const Object &other) {
  odd.clear();
  if (this != &other) {
    reset();
    import(other);
  }
  return *this;
}
Object &Object::operator<<(const Value &value) {
  if (odd.empty()) {
    odd = value.get<String>();
  } else {
    import( Object(odd, value) );
    odd.clear();
  }
  return *this;
}
Object &Object::operator<<(const Object &value) {
  import( std::string(odd),value);
  odd.clear();
  return *this;
}
size_t Object::size() const {
  return value_map_.size();
}
bool Object::empty() const {
  return value_map_.size() == 0;
}
const std::map<std::string, Value*> &Object::kv_map() const {
  return value_map_;
}
void Object::reset() {
  container::iterator i;
  for (i = value_map_.begin(); i != value_map_.end(); ++i) {
    delete i->second;
  }
  value_map_.clear();
}
bool Object::parse(std::istream &input) {
  return parse(input,*this);
}
bool Object::parse(const std::string &input) {
  std::istringstream is( input );
  return parse(is,*this);
}


Array::Array(const Array &other) {
  import(other);
}
Array::Array(const Value &value) {
  import(value);
}
void Array::import(const Array &other) {
  if (this != &other) {
    // default
    container::const_iterator
        it = other.values_.begin(),
        end = other.values_.end();
    for (/**/ ; it != end; ++it) {
      values_.push_back( new Value(**it) );
    }
  } else {
    // recursion is supported here
    import( Array(*this) );
  }
}
void Array::import(const Value &value) {
  values_.push_back( new Value(value) );
}
size_t Array::size() const {
  return values_.size();
}
bool Array::empty() const {
  return values_.size() == 0;
}
void Array::reset() {
  for (container::iterator i = values_.begin(); i != values_.end(); ++i) {
    delete *i;
  }
  values_.clear();
}
bool Array::parse(std::istream &input) {
  return parse(input,*this);
}
bool Array::parse(const std::string &input) {
  std::istringstream is(input);
  return parse(is,*this);
}
Array &Array::operator<<(const Array &other) {
  import(other);
  return *this;
}
Array &Array::operator<<(const Value &value) {
  import(value);
  return *this;
}
Array &Array::operator=(const Array &other) {
  if( this != &other ) {
    reset();
    import(other);
  }
  return *this;
}
Array &Array::operator=(const Value &value) {
  reset();
  import(value);
  return *this;
}

Value::Value(const Value &other) : type_(INVALID_) {
  import( other );
}
bool Value::empty() const {
  if( type_ == INVALID_ ) return true;
  if( type_ == STRING_ && string_value_ == 0 ) return true;
  if( type_ == ARRAY_ && array_value_ == 0 ) return true;
  if( type_ == OBJECT_ && object_value_ == 0 ) return true;
  return false;
}
bool Value::parse(std::istream &input) {
  return parse(input,*this);
}
bool Value::parse(const std::string &input) {
  std::istringstream is( input );
  return parse(is,*this);
}

}  // namespace jsonxx

using std::string;
using std::vector;
using std::pair;

#define SEQ_MAXSIZE 10000

bool operator<(pair<int, string>& p1, pair<int, string>& p2) {
  if (p1.first != p2.first) {
    return p1.first < p2.first;
  }
  return p1.second.compare(p2.second) < 0;
}

struct Seq {
  string id;
  //string v_fam;
  string v_gene;
  string v_all;
  string j_gene;
  string j_all;
  string junc;
  vector<pair<int, string>> muts;

  Seq(jsonxx::Object& o, const string& junc_query) {
    id = o.get<string>("seq_id");
    //v_fam = o.get<jsonxx::Object>("v_gene").get<string>("fam");
    v_gene = o.get<jsonxx::Object>("v_gene").get<string>("gene");
    v_all = o.get<jsonxx::Object>("v_gene").get<string>("all");
    j_gene = o.get<jsonxx::Object>("j_gene").get<string>("gene");
    j_all = o.get<jsonxx::Object>("j_gene").get<string>("all");
    junc = o.get<string>(junc_query);
    if (o.has<jsonxx::Object>("var_muts_nt")) {
      jsonxx::Array a = o.get<jsonxx::Object>("var_muts_nt").get<jsonxx::Array>("muts");
      muts.reserve(a.size());
      for (size_t i = 0; i < a.size(); i++) {
        muts.emplace_back(atoi(a.get<jsonxx::Object>(i).get<string>("loc").c_str()),
                          a.get<jsonxx::Object>(i).get<string>("mut"));
      }
    }
    sort(muts.begin(), muts.end());
  }

  Seq() {}
};

bool serialize(Seq& seq, size_t* len, char(* data)[SEQ_MAXSIZE]) {
  *len = sizeof(size_t) * 7 +
      seq.id.length() +
      seq.v_gene.length() +
      seq.v_all.length() + 
      seq.j_gene.length() +
      seq.j_all.length() +
      seq.junc.length();
  for (size_t i = 0; i < seq.muts.size(); i++) {
    *len += sizeof(int) + sizeof(size_t) + seq.muts[i].second.length();
  }
  if (*len > SEQ_MAXSIZE) {
    fprintf(stderr, "Seq too long\n");
    return false;
  }
  size_t pos = 0;
  *(size_t*)((*data) + pos) = seq.id.length();
  pos += sizeof(size_t);
  strcpy((*data) + pos, seq.id.c_str());
  pos += seq.id.length();
  *(size_t*)((*data) + pos) = seq.v_gene.length();
  pos += sizeof(size_t);
  strcpy((*data) + pos, seq.v_gene.c_str());
  pos += seq.v_gene.length();
  *(size_t*)((*data) + pos) = seq.v_all.length();
  pos += sizeof(size_t);
  strcpy((*data) + pos, seq.v_all.c_str());
  pos += seq.v_all.length();
  *(size_t*)((*data) + pos) = seq.j_gene.length();
  pos += sizeof(size_t);
  strcpy((*data) + pos, seq.j_gene.c_str());
  pos += seq.j_gene.length();
  *(size_t*)((*data) + pos) = seq.j_all.length();
  pos += sizeof(size_t);
  strcpy((*data) + pos, seq.j_all.c_str());
  pos += seq.j_all.length();
  *(size_t*)((*data) + pos) = seq.junc.length();
  pos += sizeof(size_t);
  strcpy((*data) + pos, seq.junc.c_str());
  pos += seq.junc.length();
  *(size_t*)((*data) + pos) = seq.muts.size();
  pos += sizeof(size_t);
  for (size_t i = 0; i < seq.muts.size(); i++) {
    *(int*)((*data) + pos) = seq.muts[i].first;
    pos += sizeof(int);
    *(size_t*)((*data) + pos) = seq.muts[i].second.length();
    pos += sizeof(size_t);
    strcpy((*data) + pos, seq.muts[i].second.c_str());
    pos += seq.muts[i].second.length();
  }
  if (pos != *len) {
    fprintf(stderr, "serialize error\n");
    return false;
  }
  return true;
}

void deserialize(char* data, Seq* seq) {
  size_t pos = 0;
  size_t len;
  
  len = *(size_t*)(data + pos);
  pos += sizeof(size_t);
  seq->id.clear();
  for (size_t i = 0; i < len; i++) {
    seq->id.push_back(*(data + pos++));
  }
  len = *(size_t*)(data + pos);
  pos += sizeof(size_t);
  seq->v_gene.clear();
  for (size_t i = 0; i < len; i++) {
    seq->v_gene.push_back(*(data + pos++));
  }
  len = *(size_t*)(data + pos);
  pos += sizeof(size_t);
  seq->v_all.clear();
  for (size_t i = 0; i < len; i++) {
    seq->v_all.push_back(*(data + pos++));
  }
  len = *(size_t*)(data + pos);
  pos += sizeof(size_t);
  seq->j_gene.clear();
  for (size_t i = 0; i < len; i++) {
    seq->j_gene.push_back(*(data + pos++));
  }
  len = *(size_t*)(data + pos);
  pos += sizeof(size_t);
  seq->j_all.clear();
  for (size_t i = 0; i < len; i++) {
    seq->j_all.push_back(*(data + pos++));
  }
  len = *(size_t*)(data + pos);
  pos += sizeof(size_t);
  seq->junc.clear();
  for (size_t i = 0; i < len; i++) {
    seq->junc.push_back(*(data + pos++));
  }
  len = *(size_t*)(data + pos);
  pos += sizeof(size_t);
  seq->muts.clear();
  for (size_t i = 0; i < len; i++) {
    int loc = *(int *)(data + pos);
    pos += sizeof(int);
    string mut;
    size_t len_mut = *(size_t*)(data + pos);
    pos += sizeof(size_t);
    for (size_t j = 0; j < len_mut; j++) {
      mut.push_back(*(data + pos++));
    }
    seq->muts.emplace_back(loc, mut);
  }
}


using std::string;
using std::min;
using std::max;

size_t LevenshteinDistance(string& s1, string& s2) {
  size_t dp[s1.length() + 1][s2.length() + 1];
  for (size_t i = 0; i <= s1.length(); i++) {
    dp[i][0] = i;
  }
  for (size_t i = 0; i <= s2.length(); i++) {
    dp[0][i] = i;
  }
  for (size_t i = 1; i <= s1.length(); i++) {
    for (size_t j = 1; j <= s2.length(); j++) {
      dp[i][j] = min(dp[i - 1][j] + 1,
                     dp[i][j - 1] + 1);
      dp[i][j] = min(dp[i][j],
                     dp[i - 1][j - 1] + (s1[i - 1] == s2[j - 1] ? 0 : 1));
    }
  }
  return dp[s1.length()][s2.length()];
}

double GlobalAlignment(string& s1, 
                       string& s2,
                       double match,
                       double mismatch,
                       double gap,
                       double extend) {
  assert(gap == extend);
  double dp[s1.length() + 1][s2.length() + 1];
  for (size_t i = 0; i <= s1.length(); i++) {
    dp[i][0] = gap * i;
  }
  for (size_t i = 0; i <= s2.length(); i++) {
    dp[0][i] = gap * i;
  }
  for (size_t i = 1; i <= s1.length(); i++) {
    for (size_t j = 1; j <= s2.length(); j++) {
      dp[i][j] = dp[i - 1][j - 1] + (s1[i - 1] == s2[j - 1] ? match : mismatch);
      dp[i][j] = max(dp[i][j], dp[i - 1][j] + gap);
      dp[i][j] = max(dp[i][j], dp[i][j - 1] + gap);
    }
  }
  return dp[s1.length()][s2.length()];
}

double GetLD(Seq& s1, Seq& s2) {
  if (s1.junc.length() == s2.junc.length()) {
    return s1.junc.length() - GlobalAlignment(s1.junc, s2.junc, 1, 0, -50, -50);
  } else {
    return LevenshteinDistance(s1.junc, s2.junc);
  }
}

double vCompare(Seq& s1, Seq& s2) {
  if (0 != s1.v_gene.compare(s2.v_gene)) {
    return 8;
  }
  if (0 != s1.v_all.compare(s2.v_all)) {
    return 1;
  }
  return 0;
}

double jCompare(Seq& s1, Seq& s2) {
  if (0 != s1.j_gene.compare(s2.j_gene)) {
    return 8;
  }
  if (0 != s1.j_all.compare(s2.j_all)) {
    return 1;
  }
  return 0;
}

double SharedMuts(Seq& s1, Seq& s2) {
  if (0 == s1.id.compare(s2.id)) {
    return 0.;
  }
  double bonus = 0.;
  size_t p1 = 0, p2 = 0;
  while (p1 < s1.muts.size() && p2 < s2.muts.size()) {
    if (s1.muts[p1] < s2.muts[p2]) {
      p1++;
    } else if (s2.muts[p2] < s1.muts[p1]) {
      p2++;
    } else {
      p1++;
      bonus += 0.35;
    }
  }
  return bonus;
}

double GetScore(Seq& s1, Seq& s2) {
  if (0 == s1.id.compare(s2.id)) {
    return 0.;
  }
  double LD = GetLD(s1, s2);
  double vPenalty = vCompare(s1, s2);
  double jPenalty = jCompare(s1, s2);
  double lenPenalty = fabs(0. + s1.junc.length() - s2.junc.length()) * 2;
  double editLength = (double)min(s1.junc.length(), s2.junc.length());
  double mutBonus = SharedMuts(s1, s2);
  if (mutBonus > (LD + vPenalty + jPenalty)) {
    mutBonus = (LD + vPenalty + jPenalty - 0.001);
  }
  return (LD + vPenalty + jPenalty + lenPenalty - mutBonus) / editLength;
}

#define fc_isnan(X) ((X)!=(X))


#ifndef NO_INCLUDE_FENV
#include <fenv.h>
#endif

#ifndef DBL_MANT_DIG
#error The constant DBL_MANT_DIG could not be defined.
#endif
#define T_FLOAT_MANT_DIG DBL_MANT_DIG

#ifndef LONG_MAX
#include <climits>
#endif
#ifndef LONG_MAX
#error The constant LONG_MAX could not be defined.
#endif
#ifndef INT_MAX
#error The constant INT_MAX could not be defined.
#endif

#ifndef INT32_MAX
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#endif

#ifndef HAVE_DIAGNOSTIC
#if __GNUC__ > 4 || (__GNUC__ == 4 && (__GNUC_MINOR__ >= 6))
#define HAVE_DIAGNOSTIC 1
#endif
#endif

#ifndef HAVE_VISIBILITY
#if __GNUC__ >= 4
#define HAVE_VISIBILITY 1
#endif
#endif

/* Since the public interface is given by the Python respectively R interface,
 * we do not want other symbols than the interface initalization routines to be
 * visible in the shared object file. The "visibility" switch is a GCC concept.
 * Hiding symbols keeps the relocation table small and decreases startup time.
 * See http://gcc.gnu.org/wiki/Visibility
 */
#if HAVE_VISIBILITY
#pragma GCC visibility push(hidden)
#endif

typedef int_fast32_t t_index;
#ifndef INT32_MAX
#define MAX_INDEX 0x7fffffffL
#else
#define MAX_INDEX INT32_MAX
#endif
#if (LONG_MAX < MAX_INDEX)
#error The integer format "t_index" must not have a greater range than "long int".
#endif
#if (INT_MAX > MAX_INDEX)
#error The integer format "int" must not have a greater range than "t_index".
#endif
typedef double t_float;

// self-destructing array pointer
template <typename type>
class auto_array_ptr{
private:
  type * ptr;
  auto_array_ptr(auto_array_ptr const &); // non construction-copyable
  auto_array_ptr& operator=(auto_array_ptr const &); // non copyable
public:
  auto_array_ptr()
    : ptr(NULL)
  { }
  template <typename index>
  auto_array_ptr(index const size)
    : ptr(new type[size])
  { }
  template <typename index, typename value>
  auto_array_ptr(index const size, value const val)
    : ptr(new type[size])
  {
    std::fill_n(ptr, size, val);
  }
  ~auto_array_ptr() {
    delete [] ptr; }
  void free() {
    delete [] ptr;
    ptr = NULL;
  }
  template <typename index>
  void init(index const size) {
    ptr = new type [size];
  }
  template <typename index, typename value>
  void init(index const size, value const val) {
    init(size);
    std::fill_n(ptr, size, val);
  }
  inline operator type *() const { return ptr; }
};

struct node {
  t_index node1, node2;
  t_float dist;

  /*
  inline bool operator< (const node a) const {
    return this->dist < a.dist;
  }
  */

  inline friend bool operator< (const node a, const node b) {
    return (a.dist < b.dist);
  }
};

class cluster_result {
private:
  auto_array_ptr<node> Z;
  t_index pos;

public:
  cluster_result(const t_index size)
    : Z(size)
    , pos(0)
  {}

  void append(const t_index node1, const t_index node2, const t_float dist) {
    Z[pos].node1 = node1;
    Z[pos].node2 = node2;
    Z[pos].dist  = dist;
    ++pos;
  }

  node * operator[] (const t_index idx) const { return Z + idx; }

  /* Define several methods to postprocess the distances. All these functions
     are monotone, so they do not change the sorted order of distances. */

  void sqrt() const {
    for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
      ZZ->dist = ::sqrt(ZZ->dist);
    }
  }

  void sqrt(const t_float) const { // ignore the argument
    sqrt();
  }

  void sqrtdouble(const t_float) const { // ignore the argument
    for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
      ZZ->dist = ::sqrt(2*ZZ->dist);
    }
  }

  #ifdef R_pow
  #define my_pow R_pow
  #else
  #define my_pow pow
  #endif

  void power(const t_float p) const {
    t_float const q = 1/p;
    for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
      ZZ->dist = my_pow(ZZ->dist,q);
    }
  }

  void plusone(const t_float) const { // ignore the argument
    for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
      ZZ->dist += 1;
    }
  }

  void divide(const t_float denom) const {
    for (node * ZZ=Z; ZZ!=Z+pos; ++ZZ) {
      ZZ->dist /= denom;
    }
  }
};

class doubly_linked_list {
  /*
    Class for a doubly linked list. Initially, the list is the integer range
    [0, size]. We provide a forward iterator and a method to delete an index
    from the list.

    Typical use: for (i=L.start; L<size; i=L.succ[I])
    or
    for (i=somevalue; L<size; i=L.succ[I])
  */
public:
  t_index start;
  auto_array_ptr<t_index> succ;

private:
  auto_array_ptr<t_index> pred;
  // Not necessarily private, we just do not need it in this instance.

public:
  doubly_linked_list(const t_index size)
    // Initialize to the given size.
    : start(0)
    , succ(size+1)
    , pred(size+1)
  {
    for (t_index i=0; i<size; ++i) {
      pred[i+1] = i;
      succ[i] = i+1;
    }
    // pred[0] is never accessed!
    //succ[size] is never accessed!
  }

  ~doubly_linked_list() {}

  void remove(const t_index idx) {
    // Remove an index from the list.
    if (idx==start) {
      start = succ[idx];
    }
    else {
      succ[pred[idx]] = succ[idx];
      pred[succ[idx]] = pred[idx];
    }
    succ[idx] = 0; // Mark as inactive
  }

  bool is_inactive(t_index idx) const {
    return (succ[idx]==0);
  }
};

// Indexing functions
// D is the upper triangular part of a symmetric (NxN)-matrix
// We require r_ < c_ !
#define D_(r_,c_) ( D[(static_cast<std::ptrdiff_t>(2*N-3-(r_))*(r_)>>1)+(c_)-1] )
// Z is an ((N-1)x4)-array
#define Z_(_r, _c) (Z[(_r)*4 + (_c)])

/*
  Lookup function for a union-find data structure.

  The function finds the root of idx by going iteratively through all
  parent elements until a root is found. An element i is a root if
  nodes[i] is zero. To make subsequent searches faster, the entry for
  idx and all its parents is updated with the root element.
 */
class union_find {
private:
  auto_array_ptr<t_index> parent;
  t_index nextparent;

public:
  union_find(const t_index size)
    : parent(size>0 ? 2*size-1 : 0, 0)
    , nextparent(size)
  { }

  t_index Find (t_index idx) const {
    if (parent[idx] != 0 ) { // a → b
      t_index p = idx;
      idx = parent[idx];
      if (parent[idx] != 0 ) { // a → b → c
        do {
          idx = parent[idx];
        } while (parent[idx] != 0);
        do {
          t_index tmp = parent[p];
          parent[p] = idx;
          p = tmp;
        } while (parent[p] != idx);
      }
    }
    return idx;
  }

  void Union (const t_index node1, const t_index node2) {
    parent[node1] = parent[node2] = nextparent++;
  }
};

class nan_error{};
#ifdef FE_INVALID
class fenv_error{};
#endif

/* Functions for the update of the dissimilarity array */

inline static void f_average( t_float * const b, const t_float a, const t_float s, const t_float t) {
  *b = s*a + t*(*b);
  #ifndef FE_INVALID
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wfloat-equal"
#endif
  if (fc_isnan(*b)) {
    throw(nan_error());
  }
#if HAVE_DIAGNOSTIC
#pragma GCC diagnostic pop
#endif
  #endif
}

template <typename t_members>
static void NN_chain_core(const t_index N, t_float * const D, t_members * const members, cluster_result & Z2) {
/*
    N: integer
    D: condensed distance matrix N*(N-1)/2
    Z2: output data structure

    This is the NN-chain algorithm, described on page 86 in the following book:

    Fionn Murtagh, Multidimensional Clustering Algorithms,
    Vienna, Würzburg: Physica-Verlag, 1985.
*/
  t_index i;

  auto_array_ptr<t_index> NN_chain(N);
  t_index NN_chain_tip = 0;

  t_index idx1, idx2;

  t_float size1, size2;
  doubly_linked_list active_nodes(N);

  t_float min;


  #ifdef FE_INVALID
  if (feclearexcept(FE_INVALID)) throw fenv_error();
  #endif

  for (t_index j=0; j<N-1; ++j) {
    if (NN_chain_tip <= 3) {
      NN_chain[0] = idx1 = active_nodes.start;
      NN_chain_tip = 1;

      idx2 = active_nodes.succ[idx1];
      min = D_(idx1,idx2);
      for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i]) {
        if (D_(idx1,i) < min) {
          min = D_(idx1,i);
          idx2 = i;
        }
      }
    }  // a: idx1   b: idx2
    else {
      NN_chain_tip -= 3;
      idx1 = NN_chain[NN_chain_tip-1];
      idx2 = NN_chain[NN_chain_tip];
      min = idx1<idx2 ? D_(idx1,idx2) : D_(idx2,idx1);
    }  // a: idx1   b: idx2

    do {
      NN_chain[NN_chain_tip] = idx2;

      for (i=active_nodes.start; i<idx2; i=active_nodes.succ[i]) {
        if (D_(i,idx2) < min) {
          min = D_(i,idx2);
          idx1 = i;
        }
      }
      for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i]) {
        if (D_(idx2,i) < min) {
          min = D_(idx2,i);
          idx1 = i;
        }
      }

      idx2 = idx1;
      idx1 = NN_chain[NN_chain_tip++];

    } while (idx2 != NN_chain[NN_chain_tip-2]);

    Z2.append(idx1, idx2, min);

    if (idx1>idx2) {
      t_index tmp = idx1;
      idx1 = idx2;
      idx2 = tmp;
    }

      size1 = static_cast<t_float>(members[idx1]);
      size2 = static_cast<t_float>(members[idx2]);
      members[idx2] += members[idx1];

    // Remove the smaller index from the valid indices (active_nodes).
    active_nodes.remove(idx1);

      t_float s = size1/(size1+size2);
      t_float t = size2/(size1+size2);
      for (i=active_nodes.start; i<idx1; i=active_nodes.succ[i])
        f_average(&D_(i, idx2), D_(i, idx1), s, t );
      // Update the distance matrix in the range (idx1, idx2).
      for (; i<idx2; i=active_nodes.succ[i])
        f_average(&D_(i, idx2), D_(idx1, i), s, t );
      // Update the distance matrix in the range (idx2, N).
      for (i=active_nodes.succ[idx2]; i<N; i=active_nodes.succ[i])
        f_average(&D_(idx2, i), D_(idx1, i), s, t );

  }
  #ifdef FE_INVALID
  if (fetestexcept(FE_INVALID)) throw fenv_error();
  #endif
}

#if HAVE_VISIBILITY
#pragma GCC visibility pop
#endif


#define size_(r_) ( ((r_<N) ? 1 : Z_(r_-N,3)) )

class linkage_output {
private:
  t_float * Z;

public:
  linkage_output(t_float * const Z_)
    : Z(Z_)
  {}

  void append(const t_index node1, const t_index node2, const t_float dist,
              const t_float size) {
    if (node1<node2) {
      *(Z++) = static_cast<t_float>(node1);
      *(Z++) = static_cast<t_float>(node2);
    }
    else {
      *(Z++) = static_cast<t_float>(node2);
      *(Z++) = static_cast<t_float>(node1);
    }
    *(Z++) = dist;
    *(Z++) = size;
  }
};

template <const bool sorted>
static void generate_SciPy_dendrogram(t_float * const Z, cluster_result & Z2, const t_index N) {
  // The array "nodes" is a union-find data structure for the cluster
  // identities (only needed for unsorted cluster_result input).
  union_find nodes(sorted ? 0 : N);
  if (!sorted) {
    std::stable_sort(Z2[0], Z2[N-1]);
  }

  linkage_output output(Z);
  t_index node1, node2;

  for (node const * NN=Z2[0]; NN!=Z2[N-1]; ++NN) {
    // Get two data points whose clusters are merged in step i.
    if (sorted) {
      node1 = NN->node1;
      node2 = NN->node2;
    }
    else {
      // Find the cluster identifiers for these points.
      node1 = nodes.Find(NN->node1);
      node2 = nodes.Find(NN->node2);
      // Merge the nodes in the union-find data structure by making them
      // children of a new node.
      nodes.Union(node1, node2);
    }
    output.append(node1, node2, NN->dist, size_(node1)+size_(node2));
  }
}

void linkage(const size_t N, double* matrix, double* Z) {
  //t_index* members = new t_index[N];
  auto_array_ptr<t_index> members;
  members.init(N, 1);
  cluster_result Z2(N - 1);
  NN_chain_core<t_index>(N, matrix, members, Z2);
  generate_SciPy_dendrogram<false>(Z, Z2, N);
  //delete []members;
}

#define CPY_MAX(_x, _y) ((_x > _y) ? (_x) : (_y))
#define CPY_MIN(_x, _y) ((_x < _y) ? (_x) : (_y))

#define NCHOOSE2(_n) ((_n)*(_n-1)/2)

#define CPY_BITS_PER_CHAR (sizeof(unsigned char) * 8)
#define CPY_FLAG_ARRAY_SIZE_BYTES(num_bits) (CPY_CEIL_DIV((num_bits), \
                                                          CPY_BITS_PER_CHAR))
#define CPY_GET_BIT(_xx, i) (((_xx)[(i) / CPY_BITS_PER_CHAR] >> \
                             ((CPY_BITS_PER_CHAR-1) - \
                              ((i) % CPY_BITS_PER_CHAR))) & 0x1)
#define CPY_SET_BIT(_xx, i) ((_xx)[(i) / CPY_BITS_PER_CHAR] |= \
                              ((0x1) << ((CPY_BITS_PER_CHAR-1) \
                                         -((i) % CPY_BITS_PER_CHAR))))
#define CPY_CLEAR_BIT(_xx, i) ((_xx)[(i) / CPY_BITS_PER_CHAR] &= \
                              ~((0x1) << ((CPY_BITS_PER_CHAR-1) \
                                         -((i) % CPY_BITS_PER_CHAR))))

#ifndef CPY_CEIL_DIV
#define CPY_CEIL_DIV(x, y) ((((double)x)/(double)y) == \
                            ((double)((x)/(y))) ? ((x)/(y)) : ((x)/(y) + 1))
#endif

#ifdef CPY_DEBUG
#define CPY_DEBUG_MSG(...) fprintf(stderr, __VA_ARGS__)
#else
#define CPY_DEBUG_MSG(...)
#endif


#define ISCLUSTER(_nd) ((_nd)->id >= n)
#define GETCLUSTER(_id) ((lists + _id - n))

/** The number of link stats (for the inconsistency computation) for each
    cluster. */

#define CPY_NIS 4

/** The column offsets for the different link stats for the inconsistency
    computation. */
#define CPY_INS_MEAN 0
#define CPY_INS_STD 1
#define CPY_INS_N 2
#define CPY_INS_INS 3

/** The number of linkage stats for each cluster. */
#define CPY_LIS 4

/** The column offsets for the different link stats for the linkage matrix. */
#define CPY_LIN_LEFT 0
#define CPY_LIN_RIGHT 1
#define CPY_LIN_DIST 2
#define CPY_LIN_CNT 3

void get_max_dist_for_each_cluster(const double *Z, double *max_dists, int n) {
  int *curNode;
  int ndid, lid, rid, k;
  unsigned char *lvisited, *rvisited;
  const double *Zrow;
  double max_dist;
  const int bff = CPY_FLAG_ARRAY_SIZE_BYTES(n);

  k = 0;
  curNode = (int*)malloc(n * sizeof(int));
  lvisited = (unsigned char*)malloc(bff);
  rvisited = (unsigned char*)malloc(bff);
  curNode[k] = (n * 2) - 2;
  memset(lvisited, 0, bff);
  memset(rvisited, 0, bff);
  while (k >= 0) {
    ndid = curNode[k];
    Zrow = Z + ((ndid-n) * CPY_LIS);
    lid = (int)Zrow[CPY_LIN_LEFT];
    rid = (int)Zrow[CPY_LIN_RIGHT];
    if (lid >= n && !CPY_GET_BIT(lvisited, ndid-n)) {
      CPY_SET_BIT(lvisited, ndid-n);
      curNode[k+1] = lid;
      k++;
      continue;
    }
    if (rid >= n && !CPY_GET_BIT(rvisited, ndid-n)) {
      CPY_SET_BIT(rvisited, ndid-n);
      curNode[k+1] = rid;
      k++;
      continue;
    }
    max_dist = Zrow[CPY_LIN_DIST];
    if (lid >= n) {
      max_dist = CPY_MAX(max_dist, max_dists[lid-n]);
    }
    if (rid >= n) {
      max_dist = CPY_MAX(max_dist, max_dists[rid-n]);
    }
    max_dists[ndid-n] = max_dist;
    CPY_DEBUG_MSG("i=%d maxdist[i]=%5.5f verif=%5.5f\n",
		  ndid-n, max_dist, max_dists[ndid-n]);
    k--;
  }
  free(curNode);
  free(lvisited);
  free(rvisited);
}

void form_flat_clusters_from_monotonic_criterion(const double *Z,
						 const double *mono_crit,
						 int *T, double cutoff, int n) {
  int *curNode;
  int ndid, lid, rid, k, ms, nc;
  unsigned char *lvisited, *rvisited;
  double max_crit;
  const double *Zrow;
  const int bff = CPY_FLAG_ARRAY_SIZE_BYTES(n);

  curNode = (int*)malloc(n * sizeof(int));
  lvisited = (unsigned char*)malloc(bff);
  rvisited = (unsigned char*)malloc(bff);

  /** number of clusters formed so far. */
  nc = 0;
  /** are we in part of a tree below the cutoff? .*/
  ms = -1;
  k = 0;
  curNode[k] = (n * 2) - 2;
  memset(lvisited, 0, bff);
  memset(rvisited, 0, bff);
  ms = -1;
  while (k >= 0) {
    ndid = curNode[k];
    Zrow = Z + ((ndid-n) * CPY_LIS);
    lid = (int)Zrow[CPY_LIN_LEFT];
    rid = (int)Zrow[CPY_LIN_RIGHT];
    max_crit = mono_crit[ndid-n];
    CPY_DEBUG_MSG("cutoff: %5.5f maxc: %5.5f nc: %d\n", cutoff, max_crit, nc);
    if (ms == -1 && max_crit <= cutoff) {
      CPY_DEBUG_MSG("leader: i=%d\n", ndid);
      ms = k;
      nc++;
    }
    if (lid >= n && !CPY_GET_BIT(lvisited, ndid-n)) {
      CPY_SET_BIT(lvisited, ndid-n);
      curNode[k+1] = lid;
      k++;
      continue;
    }
    if (rid >= n && !CPY_GET_BIT(rvisited, ndid-n)) {
      CPY_SET_BIT(rvisited, ndid-n);
      curNode[k+1] = rid;
      k++;
      continue;
    }
    if (ndid >= n) {
      if (lid < n) {
	if (ms == -1) {
	  nc++;
	  T[lid] = nc;
	}
	else {
	  T[lid] = nc;
	}
      }
      if (rid < n) {
	if (ms == -1) {
	  nc++;
	  T[rid] = nc;
	}
	else {
	  T[rid] = nc;
	}
      }
      if (ms == k) {
	ms = -1;
      }
    }
    k--;
  }

  free(curNode);
  free(lvisited);
  free(rvisited);  
}

void form_flat_clusters_from_dist(const double *Z, int *T,
				  double cutoff, int n) {
  double *max_dists = (double*)malloc(sizeof(double) * n);
  get_max_dist_for_each_cluster(Z, max_dists, n);
  //CPY_DEBUG_MSG("cupid: n=%d cutoff=%5.5f MD[0]=%5.5f MD[n-1]=%5.5f\n", n, cutoff, max_dists[0], max_dists[n-2]);
  form_flat_clusters_from_monotonic_criterion(Z, max_dists, T, cutoff, n);
  free(max_dists);
}

bool ParseArg(int argc, char **argv, char **inputFile, char **outputFile, size_t *n,
              bool *flagMemory); 
void linkage(const size_t N, double* matrix, double* Z);
void form_flat_clusters_from_dist(const double *Z, int *T,
				  double cutoff, int n) ;
void form_flat_clusters_maxclust_monocrit(const double *Z,
					  const double *mono_crit,
					  int *T, int n, int mc);
void get_max_dist_for_each_cluster(const double *Z, double *max_dists, int n);

struct timeval time_start;
struct timeval time_end;
struct timeval time_start_total;
struct timeval time_end_total;
#define TIC gettimeofday(&time_start, NULL);
#define TOC gettimeofday(&time_end, NULL); \
    fprintf(stderr, " %.3f sec\n", \
            (time_end.tv_sec - time_start.tv_sec) * 1. + \
            (time_end.tv_usec - time_start.tv_usec) / 1000000.);
#define TIC_TOTAL gettimeofday(&time_start_total, NULL);
#define TOC_TOTAL gettimeofday(&time_end_total, NULL); \
    fprintf(stderr, "Total Time Elapsed %.3f sec\n", \
            (time_end_total.tv_sec - time_start_total.tv_sec) * 1. + \
            (time_end_total.tv_usec - time_start_total.tv_usec) / 1000000.);

int parse_for_meminfo(char* line){
  int i = strlen(line);
  while (*line < '0' || *line > '9') line++;
  line[i-3] = '\0';
  i = atoi(line);
  return i;
}

void meminfo() {
  FILE* file = fopen("/proc/self/status", "r");
  int result = -1;
  char line[128];
  while (fgets(line, 128, file) != NULL){
    if (strncmp(line, "VmRSS:", 6) == 0){
      result = parse_for_meminfo(line);
      break;
    }
  }
  fclose(file);
  fprintf(stderr, "My Memory Usage: %d KB\n", result);
}

#define MEMINFO if (flagMemory) { meminfo(); }

using std::string;
using std::ifstream;
using std::cout;
using std::endl;
using std::vector;
using std::min;
using std::max;
using std::pair;

const double cutoff = 0.32;

bool BuildCondensedMatrix(vector<Seq>& seqs, double *scores) {
  vector<size_t> sumTable;
  sumTable.reserve(seqs.size());
  sumTable.push_back(seqs.size() - 1);
  for (size_t i = 1; i < seqs.size(); i++) {
    sumTable.push_back(sumTable.back() + seqs.size() - i - 1);
  }

  size_t nPairs = seqs.size() * (seqs.size() - 1) / 2;
  
  vector<string> ips;
  vector<int> ports;
  FILE *fSlaves = fopen("SLAVES", "r");
  size_t nSlaves = 0;

  if (fSlaves != NULL) {
    char line[1000];
    char ip[100];
    int port;
    while (fgets(line, 1000, fSlaves)) {
      if (sscanf(line, "%s%d", ip, &port) == 2 && ip[0] != '#') {
        ips.push_back(ip);
        ports.push_back(port);
        nSlaves++;
      }
    }
    fclose(fSlaves);
  }

  int sockfd[nSlaves];
  struct sockaddr_in server_addr[nSlaves];
  for (size_t i = 0; i < nSlaves; i++) {
    if ((sockfd[i] = socket(AF_INET, SOCK_STREAM, 0)) == -1) {
      fprintf(stderr, "Socket Error\n");
      return false;
    }
    memset(&server_addr[i], 0, sizeof(struct sockaddr));
    server_addr[i].sin_family = AF_INET;
    server_addr[i].sin_port = htons(ports[i]);
    struct hostent *server;
    server = gethostbyname(ips[i].c_str());
    if (server == NULL) {
      fprintf(stderr, "Find host Error\n");
      return false;
    }
    bcopy((char *)server->h_addr, (char*)&server_addr[i].sin_addr.s_addr,
          server->h_length);
    fprintf(stderr, "Connecting to %s:%d\n", ips[i].c_str(), ports[i]);
    if (connect(sockfd[i], (struct sockaddr *)&server_addr[i], 
                sizeof(struct sockaddr)) == -1) {
      fprintf(stderr, "Connect Error\n");
      return false;
    }
    fprintf(stderr, "Connected to %s:%d\n", ips[i].c_str(), ports[i]);
  }

  size_t first = 0;
  size_t second = 1;
  size_t load = nPairs / (nSlaves + 1);
  for (size_t i = 0; i < nSlaves; i++) {
    size_t nSeqs = seqs.size();
    char buf[SEQ_MAXSIZE];
    size_t len;
    
    fprintf(stderr, "Sending Assignment to %s:%d\n", ips[i].c_str(), ports[i]);
    if (send(sockfd[i], (char*)&nSeqs, sizeof(size_t), 0) == -1) {
        fprintf(stderr, "Send Error\n");
        return false;
    }
    for (size_t j = 0; j < seqs.size(); j++) {
      serialize(seqs[j], &len, &buf);
      if (send(sockfd[i], (char*)&len, sizeof(size_t), 0) == -1) {
        fprintf(stderr, "Send Error\n");
        return false;
      }
      if (send(sockfd[i], buf, len, 0) == -1) {
        fprintf(stderr, "Send Error\n");
        return false;
      }
    }
    if (send(sockfd[i], (char*)&load, sizeof(size_t), 0) == -1) {
      fprintf(stderr, "Send Error\n");
      return false;
    }
    for (size_t j = load * i; j < load * (i + 1); j++) {
      size_t p[2] = { first, second };
      if (send(sockfd[i], (char*)&p, sizeof(size_t) * 2, 0) == -1) {
        fprintf(stderr, "Send Error\n");
        return false;
      }
      second++;
      if (second == seqs.size()) {
        first++;
        second = first + 1;
      }
    }
    fprintf(stderr, "Assignment Sent to %s:%d\n", ips[i].c_str(), ports[i]);
  }
 
  fprintf(stderr, "Computing...\n");
#pragma omp parallel for 
  for (size_t i = load * nSlaves; i < nPairs; i++) {
    size_t first = std::lower_bound(sumTable.begin(), sumTable.end(), i + 1) - 
        sumTable.begin();
    size_t second = seqs.size() - 1 - (sumTable[first] - (i + 1));
    //fprintf(stderr, "%lu %lu\n", first, second);
    scores[i] =
        (float)GetScore(seqs[first], seqs[second]);
  }

  for (size_t i = 0; i < nSlaves; i++) {
    fprintf(stderr, "Receiving Result from %s:%d\n", ips[i].c_str(), ports[i]);
    for (size_t j = 0; j < load; j++) {
      if (recv(sockfd[i], (char *)(scores + load * i + j), sizeof(double), 0) != 
          sizeof(double)) {
        fprintf(stderr, "Receive Error\n");
        return false;
      }
    }
    close(sockfd[i]);
    fprintf(stderr, "Result Received from %s:%d\n", ips[i].c_str(), ports[i]);
  }
  return true;
}

bool ParseArg(int argc, char **argv, char **inputFile, char **outputFile, size_t *n,
              bool *flagMemory, char **outputFormat) {
  if (argc < 3) {
    fprintf(stderr, "argc must be no less than 3\n");
    return false;
  }
  int flag;
  *inputFile = argv[1];
  *outputFile = argv[2];
  *n = 0;
  *flagMemory = false;
  *outputFormat = NULL;
  while ((flag = getopt(argc, argv, "mn:f:")) != -1) {
    switch (flag) {
      case 'n':
        *n = atoi(optarg);
        break;
      case 'm':
        *flagMemory = true;
      case 'f':
        *outputFormat = optarg;
    }
  }
  return true;
}

int analyze(int argc, char **argv) {
  TIC_TOTAL

  char* inputFile;
  char* outputFile;
  size_t n;
  bool flagMemory;
  char* outputFormat;

  if (!ParseArg(argc, argv, &inputFile, &outputFile, &n, 
                &flagMemory, &outputFormat)) {
    return 0;
  }

  fprintf(stderr, "Loding input sequences...\n");
  TIC
  jsonxx::Array a; 
  ifstream input;
  input.open(inputFile);
  a.parse(input);
  input.close();
  vector<Seq> seqs;
  if (n == 0) {
    n = a.size();
  }
  if (a.empty()) {
    fprintf(stderr, "Error - No Sequence in the input file.\n");
    return 0;
  }
  for (size_t i = 0; i < a.size() && i < size_t(n); i++) {
    auto o = a.get<jsonxx::Object>(i);
    seqs.emplace_back(o, "junc_aa");
  }
  a.reset();
  fprintf(stderr, "Done.");
  TOC
  MEMINFO

  fprintf(stderr, "Calculating condensed distance matrix...\n");
  TIC
  double *dist_matrix = new double[seqs.size() * (seqs.size() - 1) / 2];
  if (!BuildCondensedMatrix(seqs, dist_matrix)) {
    return 0;
  }
  fprintf(stderr, "Done.");
  TOC
  MEMINFO

  fprintf(stderr, "Calculating clusters...\n");
  TIC
  double *linkage_matrix = new double[4 * (seqs.size() - 1)];
  linkage(seqs.size(), dist_matrix, linkage_matrix);
  delete []dist_matrix;
  int *flat_cluster = new int[seqs.size()];
  form_flat_clusters_from_dist(linkage_matrix, flat_cluster,
				  cutoff, seqs.size());
  delete []linkage_matrix;
  fprintf(stderr, "Done.");
  TOC
  MEMINFO

  fprintf(stderr, "Outputing clusters...\n");
  TIC
  FILE *output = fopen(outputFile, "w");
  if (outputFormat != NULL && strcmp(outputFormat, "seqs") == 0) {
    vector<vector<size_t>> dict;
    dict.resize(seqs.size());
    for (size_t i = 0; i < seqs.size(); i++) {
      dict[flat_cluster[i]].push_back(i);
    }
    for (size_t i = 0; i < seqs.size(); i++) {
      if (dict[i].size() < 2) continue;
      fprintf(output, "#lineage_v0_%lu\n", i);
      for (size_t j = 0; j < dict[i].size(); j++) {
        fprintf(output, ">%s\n%s\n", seqs[dict[i][j]].id.c_str(),
                seqs[dict[i][j]].junc.c_str());
      }
      fprintf(output, "\n");
    }
  } else {
    for (size_t i = 0; i < seqs.size(); i++) {
      fprintf(output, "%d\n", flat_cluster[i]);
    }
  }
  fclose(output);
  delete []flat_cluster;
  fprintf(stderr, "Done.");
  TOC
  MEMINFO

  TOC_TOTAL
  
  return 0;
}

#define BUFSIZE 20000
#define SERVERPORT 8888
#define MAXCONN_NUM 1

struct Assignment {
  Seq* seqs;
  pair<size_t, size_t>* pairs;
  size_t n_seqs;
  size_t n_pairs;
  double *scores;
  Assignment() : seqs(NULL), pairs(NULL), n_seqs(0), n_pairs(0), scores(NULL) {}
  ~Assignment() {
    if (seqs != NULL) {
      delete []seqs;
    }
    if (pairs != NULL) {
      delete []pairs;
    }
    if (scores != NULL) {
      delete []scores;
    }
  }
};

int slave(int argc, char **argv) {
  int port = SERVERPORT;
  if (argc >= 2) {
    port = atoi(argv[1]);
  }
  char buf[BUFSIZE];
  int sockfd, new_fd, numbytes;
  struct sockaddr_in server_addr;
  struct sockaddr_in client_addr;
  socklen_t sin_size;
  if ((sockfd = socket(AF_INET, SOCK_STREAM, 0)) == -1) {
    fprintf(stderr, "Socket Error\n");
    return 0;
  }
  memset(&client_addr, 0, sizeof(struct sockaddr));
  server_addr.sin_family = AF_INET;
  server_addr.sin_port = htons(port);
  server_addr.sin_addr.s_addr = INADDR_ANY;
  if (bind(sockfd, (struct sockaddr*)& server_addr, sizeof(struct sockaddr)) ==
      -1) {
    fprintf(stderr, "Bind Error\n");
    return 0;
  }
  if (listen(sockfd, MAXCONN_NUM) == -1) {
    fprintf(stderr, "Listen Error\n");
    return 0;
  }
  while (1) {
    fprintf(stderr, "Ready for New Assignment\n");
    sin_size = sizeof(struct sockaddr_in);
    if ((new_fd = accept(sockfd, (struct sockaddr*)&client_addr, &sin_size)) ==
        -1) {
      fprintf(stderr, "Accept Error\n");
      continue;
    }
    fprintf(stderr, "New Connection\n");
    
    Assignment assignment;
    
    if ((numbytes = recv(new_fd, buf, sizeof(size_t), 0)) != sizeof(size_t)) {
      fprintf(stderr, "Receive Data Error\n");
      return 0;
    }
    assignment.n_seqs = *(size_t*)buf;
    assignment.seqs = new Seq[assignment.n_seqs];
    for (size_t i = 0; i < assignment.n_seqs; i++) {
      if ((numbytes = recv(new_fd, buf, sizeof(size_t), 0)) != sizeof(size_t)) {
        fprintf(stderr, "Receive Data Error\n");
        return 0;
      }
      size_t len = *(size_t*)buf;
      if (len > BUFSIZE) {
        fprintf(stderr, "Seq too long\n");
        return 0;
      }
      if ((numbytes = recv(new_fd, buf, len, 0)) != len) {
        fprintf(stderr, "Receive Data Error\n");
        return 0;
      }
      deserialize(buf, &assignment.seqs[i]);
    }
    if ((numbytes = recv(new_fd, buf, sizeof(size_t), 0)) != sizeof(size_t)) {
      fprintf(stderr, "Receive Data Error\n");
      return 0;
    }
    assignment.n_pairs = *(size_t*)buf;
    assignment.pairs = new pair<size_t, size_t>[assignment.n_pairs];
    for (size_t i = 0; i < assignment.n_pairs; i++) {
      if ((numbytes = recv(new_fd, buf, sizeof(size_t) * 2, 0)) !=
          sizeof(size_t) * 2) {
        fprintf(stderr, "Receive Data Error\n");
        return 0;
      }
      assignment.pairs[i].first = *(size_t*)buf;
      assignment.pairs[i].second = *(((size_t*)buf) + 1);
    }

    assignment.scores = new double[assignment.n_pairs];
    fprintf(stderr, "Assignment Received. Start Computing.\n");
#pragma omp parallel for 
    for (size_t i = 0; i < assignment.n_pairs; i++) {
      assignment.scores[i] =
        (float)GetScore(assignment.seqs[assignment.pairs[i].first],
                        assignment.seqs[assignment.pairs[i].second]);
    }
    fprintf(stderr, "Assignment Finished. Send back to Master.\n");
    if (send(new_fd, (char *)(assignment.scores), assignment.n_pairs *
             sizeof(double), 0) == -1) {
      fprintf(stderr, "Send Error\n");
      return 0;
    } 
    close(new_fd);
  }
  close(sockfd);
  return 0;
}

int main(int argc, char** argv) {
#ifdef SLAVE
  return slave(argc, argv);
#else
  return analyze(argc, argv);
#endif 
}

